#### Set-Up ####
rm(list = ls())
Sys.info()[7]

# condition_data <- 0
# use_full_sample <- 0
expand <- 1
K <- 1000 # number of times to expand
N_bootstrap <- 1000
pa <- 0.005
pfp <- 0.01
pfn <- 0.01
J <- 1 

# Set file paths 
DATA <- "../external_links/real"
TEMP <- "../temp"
RESULTS <- "../output"
MODEL_RESULTS <- "../external/estimate_w_xs" # for separate dir

# Go to code folder
setwd("./")
sink("bootstrap_cf.txt")
folder <- paste0(TEMP, "/cf_bs_iters")
if (file.exists(folder)) {
  cat("The folder already exists")
} else {
  dir.create(folder)
}
Sys.info()[7]
Sys.time()

# Packages
library(MASS)
library(tmvtnorm)
library(haven)
library(dplyr)
library(tidyr)
library(nloptr)
library(corpcor)
library(stargazer)
library(pracma)
library(foreign)
library(writexl)
library(broom)
library(png)
library(tidyverse)

###################################################################################################################################################################
#### Load Data + Assumptions ####
analysis <- read_dta(paste0(DATA, "/model_sample.dta")) %>%
  mutate(oop_i = ifelse(is.na(oop_i)==TRUE, Inf, oop_i)) %>%
  mutate(age_25_34 = ifelse(age >= 25 & age < 35, 1, 0)) %>%
  mutate(age_35_plus = ifelse(age>35,1,0)) %>%
  mutate(some_college = ifelse(educ == 1 | is.na(educ), 0, 1)) %>%
  mutate(missing_college = ifelse(is.na(educ), 1, 0)) %>%
  mutate(any_prev_birth_issue = ifelse(prev_concern == 1 | prev_any_q_icd == 1 | prev_mis_still == 1, 1,0)) %>%
  mutate(inc_quartile_4 = ifelse(inc_quartile != 4 | is.na(inc_quartile), 0, 1)) %>%
  mutate(full_college = ifelse(educ != 3 | is.na(educ), 0, 1))

print("Begin time")
start_time <- Sys.time()
print(start_time)

#### Load Necessary Model Functions ####
source("r_functions/production_model_xs_new.R")
source("r_functions/counterfactual_functions.R")

# Function that calculates decisions for each observation i 
  # NEW: WITHOUT DRAWING A+C (i.e., taking a+c as given)
  # NEW: TAKING REC AS GIVEN 
  # NEW: ADD POSSIBLE OOP cost FOR INVASIVE = oop2_i 
  # NEW: allow mother's prior to be different from the truth (where truth=p_i=KUB)

#### Simulate data based on parameter estimates + model ####
analysis <- analysis %>%
  mutate(full_college_norm = (full_college - mean(full_college))) %>%
  mutate(fc_n = (full_college - mean(full_college))) %>%
  mutate(inc_quartile_4_norm = (inc_quartile_4 - mean(inc_quartile_4))) %>%
  mutate(married_norm = (married - mean(married))) %>%
  mutate(mom_foreign_norm = (mom_foreign - mean(mom_foreign))) %>%
  mutate(any_prev_birth_issue_norm = (any_prev_birth_issue - mean(any_prev_birth_issue))) %>%
  mutate(prev_kids_norm = (dv_prev_kids - mean(dv_prev_kids))) %>%
  mutate(age_35_plus_norm = (age_35_plus - mean(age_35_plus)))
f_a = ~ 0 + full_college_norm + inc_quartile_4_norm + married_norm + mom_foreign_norm + any_prev_birth_issue_norm + age_35_plus_norm + prev_kids_norm
f_c = ~ 0 + full_college_norm + inc_quartile_4_norm + married_norm + mom_foreign_norm + any_prev_birth_issue_norm + age_35_plus_norm + prev_kids_norm
x_vars <- c("full_college_norm", "inc_quartile_4_norm", "married_norm", "mom_foreign_norm", "any_prev_birth_issue_norm", "age_35_plus_norm", "prev_kids_norm")

## calculate true positive indicator for chrom. ab. = live birth with chrom. ab. (using actual, not simulated, data), and simulate positive indicator for counterfactuals

actual <- analysis

bootstrap_iters_cf0 <- matrix(NA, nrow = N_bootstrap, ncol = 18)
bootstrap_iters_cf1 <- matrix(NA, nrow = N_bootstrap, ncol = 18)
bootstrap_iters_cf2 <- matrix(NA, nrow = N_bootstrap, ncol = 18)
bootstrap_iters_cf4 <- matrix(NA, nrow = N_bootstrap, ncol = 18)
bootstrap_iters_cf5 <- matrix(NA, nrow = N_bootstrap, ncol = 18)
bootstrap_iters_cf7 <- matrix(NA, nrow = N_bootstrap, ncol = 18)
bootstrap_iters_cf8 <- matrix(NA, nrow = N_bootstrap, ncol = 18)
bootstrap_iters_cf9 <- matrix(NA, nrow = N_bootstrap, ncol = 18)

if (file.exists(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf0.csv"))) {
  file.remove(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf0.csv"))
  file.remove(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf1.csv"))
  file.remove(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf2.csv"))
  file.remove(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf4.csv"))
  file.remove(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf5.csv"))
  file.remove(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf7.csv"))
  file.remove(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf8.csv"))
  file.remove(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf9.csv"))
}

for (x in 0:N_bootstrap) {

  # set bootstrap sample
  if (x > 0) {
    # bootstrap_index <- sample(1:nrow(actual), size = smp_size, replace = TRUE)
    # bootstrap_sample <- actual[bootstrap_index,]
    actual <- read.csv(file=paste0(TEMP, "/bs_samples/bs_sample_", x, ".csv"), na=".")
  }

  # get estimated parameters for this sample
  if (x > 0) {
      est_par <- read.csv(paste0(RESULTS, "/bootstrap_se.csv")) %>%
        as.matrix()
      est_par <- est_par[x:x,2:21]
    } else {
      est_par <- read.csv(paste0(MODEL_RESULTS, "/estimate_w_xs/results_no_xs.csv")) %>%
        as.matrix() 
      est_par <- est_par[-21:-22,]
    }
  print(est_par)

  #Step 0a: expand data K times
  if (expand == 1) {
    actual_1 <- actual
    if (K > 1) {
      for (i in 1:(K-1)) {
        actual <- rbind(actual, actual_1)
      }
    }
    rm(i, actual_1)
  }

  print(paste0("bootstrap iteration ", x))
  seed <- x + 321
  set.seed(seed)
  smp_size <- as.numeric(nrow(analysis))

 

  ###################################################################################################################################################################
  #Step 0b: simulate positive indicator for chromAb 
  set.seed(1292837)
  actual <- actual %>% 
    mutate(draw = runif(nrow(.), min=0, max=1)) %>%
    mutate(sim_positive = ifelse(draw <= p_i, 1, 0)) %>% 
    select(-draw)

  #Step 0c: simulate result of NIPT test, if taken (using simulated chromAb indicator)
  # added new mutate(sim_nipt_result) before select(-draw) to address dropping pregnancies in conditional counterfactual
  actual <- actual %>% 
    mutate(draw = runif(nrow(.), min=0, max=1)) %>%
    mutate(sim_nipt_result = ifelse(sim_positive == 1,
                                 ifelse(draw > pfn, 1, 0),
                                 ifelse(draw <= pfp, 1, 0))) %>%
  select(-draw)

  ###### CHECKS FOR SIM_NIPT_RESULT ######
  print("sim_nipt_results after changing procedure")
  print("table of actual$sim_positive")
  table(actual$sim_positive)
  pos_ca <- actual %>%
    filter(sim_positive == 1)
  neg_test_ca <- pos_ca %>%
    filter(sim_nipt_result == 0)
  print("pos_ca$sim_nipt_result")
  table(pos_ca$sim_nipt_result)
  pos_ca_count <- nrow(pos_ca)
  neg_test_pos_ca <- nrow(neg_test_ca)
  false_neg_rate <- neg_test_pos_ca/pos_ca_count
  print("false negative rate")
  print(false_neg_rate)
  
  neg_ca <- actual %>%
    filter(sim_positive == 0)
  pos_test_no_ca <- neg_ca %>%
    filter(sim_nipt_result == 1)
  print("neg_ca$sim_nipt_result")
  table(neg_ca$sim_nipt_result)
  neg_ca_count <- nrow(neg_ca)
  pos_test_neg_ca <- nrow(pos_test_no_ca)
  false_pos_rate <- pos_test_neg_ca/neg_ca_count
  print("false positive rate")
  print(false_pos_rate)
  
  rm(pos_ca, neg_test_ca, neg_ca, pos_test_no_ca)
  table(actual$sim_nipt_result, useNA = "ifany")

  # Step 1: simulate K (a_i, c_i) draws per pregnancy (using parameter estimates) + calculate predicted testing decisions
  # note: condition on the data = discard (a_i, c_i) draws that are rejected by the observed testing decisions
  draw <- draw_a_c_counter_xs(par=est_par, data=actual, f_a=f_a, f_c=f_c, x_vars = x_vars)
  draw <- as.data.frame(draw) %>%
    mutate(oop2_i = 0)
  simulated <- model_decisions_counterfactuals(par=est_par, data=draw, prior=quo(p_i), setup = 1) %>%
    as.data.frame() %>%
    dplyr::select(a_i, c_i, pred_nipt, pred_invasive)
  all <- cbind(actual, simulated)
  all <- all %>% mutate(sim_wgt = 1) 

  # Step 2: for each observation, randomly draw whether invasive testing will cause a miscarriage (unknown to mother)
  all <- all %>%
    mutate(draw = runif(nrow(.), min=0, max=1)) %>%
    mutate(sim_inv_mis = ifelse(draw <= pa, 1, 0)) %>% 
    select(-draw) 
    
  # Step 3: clean up
  # only keep necessary variables 
  # + create new oop cost for invasive variable (=0 for all, which is actual)
  # + create new rec variable (= actual, for now)
  all <- all %>%
    dplyr::select(pregnancy, did_nipt, did_invasive, live_birth, p_i, wave, oop_i, age, ave_kub_age, sim_positive, a_i, c_i, sim_wgt, sim_inv_mis, sim_nipt_result, fetus_risk) %>%
    mutate(oop2_i = 0) %>%
    mutate(rec_i = ifelse(fetus_risk <= 200, 1, 0)) %>% 
    mutate(low = ifelse(fetus_risk > 200, 1, 0)) %>%
    mutate(med = ifelse(fetus_risk <= 200 & p_i >= 51, 1, 0)) %>%
    mutate(high = ifelse(fetus_risk <= 51, 1, 0))

  print("Weighted share of accepted draws with c_i greater than a_i" )
  all <- all %>% mutate(c_g_a = ifelse(c_i > a_i, 1, 0))
  print(weighted.mean(all$c_g_a, all$sim_wgt))

  ###################################################################################################################################################################
  #### Calc counterfactual decisions and outcomes ####
    
  # n pregnancies
  n_preg <- nrow(analysis)
  share_low <- weighted.mean(all$low, all$sim_wgt)
  share_med <- weighted.mean(all$med, all$sim_wgt)
  share_high <- weighted.mean(all$high, all$sim_wgt)
  sum_wgts <- sum(all$sim_wgt)
  
  # Set cost to government of tests (in dollars)
  g_cost_inv <- 1248.50 # 11,000 SEK 
  g_cost_nipt <- 567.5 # 5000 SEK
  g_cost_nt <- 174 #1,500 SEK
    
  # Counterfactual 0 = model predictions
  # Implement settings
  set <- all
  # Predict decisions
  dec <- model_decisions_counterfactuals(par=est_par, data=set, prior=quo(p_i))
  # Calculate outcomes 
  out <- outcomes(par=est_par, data=dec, prior=quo(p_i))
  # Calculate shares of outcomes 
  cf0 <- outcomes_shares(data=out)
  if (x > 0) {
    # bootstrap_iters_cf0[x,] <- cf0
    cf <- as.matrix(t(cf0))
    write.table(cf, file = paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf0.csv"), sep=",", append=TRUE, row.names=FALSE, col.names=FALSE)
  }

  rm(set, dec, out) 
  
  ######################################################################
  # Counterfactual 1 = mothers already know whether fetus has chrom ab
  # settings: SPECIAL - don't use "outcomes" function to calc outcomes
  #   no testing
  #   abort if positive==1 & a_i > c_i
    
  # Predict decisions
  dec <- all %>% 
    mutate(pred_nipt = 0) %>%
    mutate(pred_invasive = 0)
  
  # Calculate outcomes 
  out <- dec %>%
    mutate(sim_abort = ifelse((sim_positive==1 & (a_i > c_i)), 1, 0)) %>% 
    mutate(sim_live = ifelse(sim_abort==1, 0, 1)) %>% 
    mutate(sim_mis = 0) %>% 
    mutate(sim_terminate = sim_abort + sim_mis) %>%
    mutate(sim_nipt = pred_nipt) %>% # just renaming 
    mutate(sim_invasive = pred_invasive) %>% # just renaming
    mutate(sim_util = ifelse(sim_live == 0, a_i,
                             ifelse(sim_positive == 1, c_i, 0))) %>% 
    mutate(sim_util = ifelse(pred_nipt == 0, sim_util - (oop2_i * pred_invasive), sim_util - oop_i - oop2_i * (pred_invasive))) %>%
    dplyr::select(sim_wgt, sim_positive, sim_nipt, sim_invasive, sim_live, sim_abort, sim_mis, sim_terminate, sim_util, a_i, c_i, oop_i, oop2_i, p_i, fetus_risk)  
  
  out_low <- out %>% filter(fetus_risk >200)
  out_med <- out %>% filter(fetus_risk <= 200 & fetus_risk > 50)
  out_high <- out %>% filter(fetus_risk <= 50)
  
  # Calculate shares of outcomes 
  cf1 <- outcomes_shares(data=out)
  cf1_low <- outcomes_shares_split(data=out_low)
  cf1_med <- outcomes_shares_split(data=out_med)
  cf1_high <- outcomes_shares_split(data=out_high)
  
  rm(dec, out, out_low, out_med, out_high)
  
  # Calculate cost to govt from testing
  cf1[14] = 0
  cf1[15] = 0
  cf1[16] = 0
  if (x > 0) {
    cf <- as.matrix(t(cf1))
    write.table(cf, file = paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf1.csv"), sep=",", append=TRUE, row.names=FALSE, col.names=FALSE)
  }

  ######################################################################
  # Counterfactual 2 = everyone must not get any testing
  # settings: 
  #   oop_i = Inf (for all)
  #   oop2_i = Inf (for all) 
  #   rec_i = 0 (for all)
  #   prior = p_i (i.e., prior = truth = KUB score)
    
  # Implement settings
  set <- all %>%
    mutate(oop_i = Inf) %>% 
    mutate(oop2_i = Inf)
  
  # Predict decisions
  dec <- model_decisions_counterfactuals(par=est_par, data=set, prior=quo(p_i))
  
  # Calculate outcomes 
  out <- outcomes(par=est_par, data=dec, prior=quo(p_i))
  out_low <- out %>% filter(fetus_risk >200)
  out_med <- out %>% filter(fetus_risk <= 200 & fetus_risk > 50)
  out_high <- out %>% filter(fetus_risk <= 50)
  
  # Calculate shares of outcomes 
  cf2 <- outcomes_shares(data=out)
  cf2_low <- outcomes_shares_split(data=out_low)
  cf2_med <- outcomes_shares_split(data=out_med)
  cf2_high <- outcomes_shares_split(data=out_high)
  
  rm(set, dec, out, out_low, out_med, out_high)
  
  # normalize consumer surplus to be relative to no testing CS
  surplus_norm <- cf2[17]
  cf1[17] = cf1[17] - surplus_norm
  cf2[17] = cf2[17] - surplus_norm
  if (x > 0) {
    cf <- as.matrix(t(cf2))
    write.table(cf, file = paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf2.csv"), sep=",", append=TRUE, row.names=FALSE, col.names=FALSE)
  }
    
  ######################################################################
  # Counterfactual 4 = everyone must not get NIPT; invasive offered to everyone (and free to them)
  # settings: 
  #   oop_i = Inf (for all)
  #   oop2_i = 0 (for all)
  #   rec_i = 0 (for all)
  #   prior = p_i (i.e., prior = truth = KUB score)
    
  # Implement settings
  set <- all %>%
    mutate(oop_i = Inf) %>% 
    mutate(oop2_i = 0)  

  # Predict decisions
  dec <- model_decisions_counterfactuals(par=est_par, data=set, prior=quo(p_i))
  
  # Calculate outcomes 
  out <- outcomes(par=est_par, data=dec, prior=quo(p_i))
  out_low <- out %>% filter(fetus_risk > 200)
  out_med <- out %>% filter(fetus_risk <= 200 & fetus_risk > 50)
  out_high <- out %>% filter(fetus_risk <= 50)
  
  # Calculate shares of outcomes 
  cf4 <- outcomes_shares(data=out)
  cf4_low <- outcomes_shares_split(data=out_low)
  cf4_med <- outcomes_shares_split(data=out_med)
  cf4_high <- outcomes_shares_split(data=out_high)
  
  # save decisions for plots of testing rates bycounterfactual
  write.dta(out, file=paste0(TEMP, "/inv_only_cf_out.dta"))
  rm(set, dec, out, out_low, out_med, out_high)
  
  # normalize consumer surplus to be relative to no testing CS
  cf4[17] = cf4[17] - surplus_norm
  if (x > 0) {
    cf <- as.matrix(t(cf4))
    write.table(cf, file = paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf4.csv"), sep=",", append=TRUE, row.names=FALSE, col.names=FALSE)
  }
    
  ######################################################################
  # Counterfactual 5 = everyone gets free testing
  # settings: 
  #   oop_i = 0 (for all)
  #   oop2_i = 0 (for all)
  #   rec_i = 0 (for all)
  #   prior = p_i (i.e., prior = truth = KUB score)
    
  # Implement settings
  set <- all %>%
    mutate(oop_i = 0) %>% 
    mutate(oop2_i = 0)
    
  # Predict decisions
  dec <- model_decisions_counterfactuals(par=est_par, data=set, prior=quo(p_i))

  # Calculate outcomes 
  out <- outcomes(par=est_par, data=dec, prior=quo(p_i))
  out_low <- out %>% filter(fetus_risk > 200)
  out_med <- out %>% filter(fetus_risk <= 200 & fetus_risk > 50)
  out_high <- out %>% filter(fetus_risk <= 50)
  
  # Calculate shares of outcomes 
  cf5 <- outcomes_shares(data=out)
  cf5_low <- outcomes_shares_split(data=out_low)
  cf5_med <- outcomes_shares_split(data=out_med)
  cf5_high <- outcomes_shares_split(data=out_high)
  
  # save decisions for plots of testing rates by counterfactual
  write.dta(out, file=paste0(TEMP, "/free_nipt_cf_out.dta"))
 
  rm(set, dec, out, out_low, out_med, out_high)  
  
  # normalize consumer surplus to be relative to no testing CS
  cf5[17] = cf5[17] - surplus_norm
  if (x > 0) {
    cf <- as.matrix(t(cf5))
    write.table(cf, file = paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf5.csv"), sep=",", append=TRUE, row.names=FALSE, col.names=FALSE)
  }
    
  ######################################################################
  # Counterfactual 7 = everyone pay for own NIPT but get free invasive
  # settings: 
  #   oop_i = 5000 SEK ~= 567.50 USD (for all)
  #   oop2_i = 0 (for all)
  #   rec_i = 0 (for all)
  #   prior = p_i (i.e., prior = truth = KUB score)
    
  # Implement settings
  set <- all %>%
    mutate(oop_i = 567.50) %>% 
    mutate(oop2_i = 0)
  
  # Predict decisions
  dec <- model_decisions_counterfactuals(par=est_par, data=set, prior=quo(p_i))
  
  # Calculate outcomes 
  out <- outcomes(par=est_par, data=dec, prior=quo(p_i))
  out_low <- out %>% filter(fetus_risk > 200)
  out_med <- out %>% filter(fetus_risk <= 200 & fetus_risk > 50)
  out_high <- out %>% filter(fetus_risk <= 50)
  
  # Calculate shares of outcomes 
  cf7 <- outcomes_shares(data=out)
  cf7_low <- outcomes_shares_split(data=out_low)
  cf7_med <- outcomes_shares_split(data=out_med)
  cf7_high <- outcomes_shares_split(data=out_high)
  
  # save decisions for plots of testing rates by
  write.dta(out, file=paste0(TEMP, "/oop_cf_out.dta"))

  rm(set, dec, out, out_low, out_med, out_high)
  
  # normalize consumer surplus to be relative to no testing CS
  cf7[17] = cf7[17] - surplus_norm

  if (x > 0) {
    cf <- as.matrix(t(cf7))
    write.table(cf, file = paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf7.csv"), sep=",", append=TRUE, row.names=FALSE, col.names=FALSE)
  }
    
  ######################################################################
  # Counterfactual 8 = free invasive for all, free NIPT for 1/200 <= KUB <= 1/51, pay for own NIPT otherwise
  # settings: 
  #   oop_i = 0 for 1/200 <= KUB <= 1/51, oop_i = 5000 SEK ~= 567.50 USD for others
  #   oop2_i = 0 (for all)
  #   rec_i = 0 (for all)
  #   prior = p_i (i.e., prior = truth = KUB score)
    
  # Implement settings
  set <- all %>%
    mutate(oop_i = ifelse((fetus_risk <= 200) & (fetus_risk > 50), 0, 567.50)) %>% 
    mutate(oop2_i = 0)
  
  # Predict decisions
  dec <- model_decisions_counterfactuals(par=est_par, data=set, prior=quo(p_i))
  
  # Calculate outcomes 
  out <- outcomes(par=est_par, data=dec, prior=quo(p_i))
  
  # Calculate shares of outcomes 
  cf8 <- outcomes_shares(data=out)
  rm(set, dec, out)
  
  # normalize consumer surplus to be relative to no testing CS
  cf8[17] = cf8[17] - surplus_norm

  if (x > 0) {
    cf <- as.matrix(t(cf8))
    write.table(cf, file = paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf8.csv"), sep=",", append=TRUE, row.names=FALSE, col.names=FALSE)
  }
  
  ######################################################################
  # Counterfactual 9 = free invasive for all, free NIPT for 1/200 <= KUB, pay for own NIPT otherwise
  # settings: 
  #   oop_i = 0 for 1/200 <= KUB, oop_i = 5000 SEK ~= 567.50 USD for others
  #   oop2_i = 0 (for all)
  #   rec_i = 0 (for all)
  #   prior = p_i (i.e., prior = truth = KUB score)
    
  # Implement settings
  set <- all %>%
    mutate(oop_i = ifelse(fetus_risk <= 200, 0, 567.50)) %>% 
    mutate(oop2_i = 0) 
  
  # Predict decisions
  dec <- model_decisions_counterfactuals(par=est_par, data=set, prior=quo(p_i))
  
  # Calculate outcomes 
  out <- outcomes(par=est_par, data=dec, prior=quo(p_i))
  
  # Calculate shares of outcomes 
  cf9 <- outcomes_shares(data=out)
  
  # save decisions for plots of testing rates by 
  write.dta(out, file=paste0(TEMP, "/g_200_cf_out.dta"))
  
  rm(set, dec, out)
  
  # normalize consumer surplus to be relative to no testing CS
  cf9[17] = cf9[17] - surplus_norm
  if (x > 0) {
    cf <- as.matrix(t(cf9))
    write.table(cf, file = paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf9.csv"), sep=",", append=TRUE, row.names=FALSE, col.names=FALSE)
  }

  ######################################################################
  # Full counterfactuals table
  if (x == 0) {
    out_cf <- rbind(cf0, cf1, cf2, cf4, cf5, cf7, cf8, cf9)
    rownames(out_cf) <- c("Model", "First-best", "No testing", "Invasive only", "Free NIPT", "Out-of-pocket NIPT", "[1:51-1:200]", " > 1:201")
    colnames(out_cf) <- c(c("Any testing", "NIPT only", "Invasive only", "NIPT followed by invasive", 
                            "Live birth", "live, chrom ab", "live, chrom ab, a_i > c_i", "live, chrom ab, a_i < c_i", "live, no chrom ab", 
                            "Terminated", "term, chrom ab, a_i > c_i", "term, no chrom ab", 
                            "Bad outcome", 
                            "G cost (p.c.)", "G cost, NIPT", "G cost, Inv", "Consumer surplus", "Share NIPT oop"))
    print(t(out_cf))
    stargazer(t(out_cf), digits = 2)
    write.csv(out_cf, file=paste0(RESULTS, "/cf_outcomes.csv"), na=".", row.names=TRUE)
  }
}

# Full counterfactual SE table
bootstrap_iters_cf0 <- read.csv(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf0.csv"), header=FALSE)
bootstrap_iters_cf1 <- read.csv(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf1.csv"), header=FALSE)
bootstrap_iters_cf2 <- read.csv(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf2.csv"), header=FALSE)
bootstrap_iters_cf4 <- read.csv(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf4.csv"), header=FALSE)
bootstrap_iters_cf5 <- read.csv(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf5.csv"), header=FALSE)
bootstrap_iters_cf7 <- read.csv(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf7.csv"), header=FALSE)
bootstrap_iters_cf8 <- read.csv(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf8.csv"), header=FALSE)
bootstrap_iters_cf9 <- read.csv(paste0(TEMP, "/cf_bs_iters/bootstrap_iters_cf9.csv"), header=FALSE)
sd0 <- apply(bootstrap_iters_cf0, 2, sd)
sd1 <- apply(bootstrap_iters_cf1, 2, sd)
sd2 <- apply(bootstrap_iters_cf2, 2, sd)
sd4 <- apply(bootstrap_iters_cf4, 2, sd)
sd5 <- apply(bootstrap_iters_cf5, 2, sd)
sd7 <- apply(bootstrap_iters_cf7, 2, sd)
sd8 <- apply(bootstrap_iters_cf8, 2, sd)
sd9 <- apply(bootstrap_iters_cf9, 2, sd)
out_sd <- rbind(sd0, sd1, sd2, sd4, sd5, sd7, sd8, sd9)
rownames(out_sd) <- c("Model", "First-best", "No testing", "Invasive only", "Free NIPT", "Out-of-pocket NIPT", "[1:51-1:200]", " > 1:201")
colnames(out_sd) <- c(c("Any testing", "NIPT only", "Invasive only", "NIPT followed by invasive", 
                        "Live birth", "live, chrom ab", "live, chrom ab, a_i > c_i", "live, chrom ab, a_i < c_i", "live, no chrom ab", 
                        "Terminated", "term, chrom ab, a_i > c_i", "term, no chrom ab", 
                        "Bad outcome", 
                        "G cost (p.c.)", "G cost, NIPT", "G cost, Inv", "Consumer surplus", "Share NIPT oop"))
print(paste0("number of bootstrap iterations: ", N_bootstrap))
print(t(out_sd))
stargazer(t(out_sd), digits = 2)
write.csv(out_sd, file=paste0(RESULTS, "/cf_outcomes_se.csv"), na=".", row.names=TRUE)

sink()
quit()

###################################################################################################################################################################
### Counterfactuals of varying lower bound for NIPT coverage ####
N_grid <- 50
grid <- seq(from = 0, to = 1000, by = (1000 - 0) / N_grid)
vary_lb_cf <- matrix(NA, nrow = N_grid + 1, ncol = 10)

# loop over values of lower bound and save bad outcomes/govt_cost
for (i in 1:(N_grid + 1)){
  print(i)
  lb <- grid[i]
  if (grid[i] == 0) {
    lb <- 2
  }
  vary_lb_cf[i,1] = lb
  # Implement settings
  set <- all %>%
    mutate(oop_i = ifelse(fetus_risk <= lb, 0, 567.50)) %>% 
    mutate(oop2_i = 0) 
  if ((lb) == 2) {
    set <- all %>% mutate(oop_i =  567.50)
  }
  if ((lb) == 1000) {
    set <- all %>% mutate(oop_i =  0)
  }
  
  # Predict decisions
  dec <- model_decisions_counterfactuals(par=est_par, data=set, prior=quo(p_i))
  
  # Calculate outcomes 
  out <- outcomes(par=est_par, data=dec, prior=quo(p_i))
  
  # Calculate shares of outcomes 
  out_shares <- outcomes_shares(data=out)
  
  # save bad outcome shares (+ separately by type of bad outcome), gcost, and testing
  vary_lb_cf[i,2] = out_shares[13] # bad outcome share
  vary_lb_cf[i,3] = out_shares[14] # gcost
  vary_lb_cf[i,4] = out_shares[1] # subt
  vary_lb_cf[i,5] = out_shares[2] # nipt only
  vary_lb_cf[i,6] = out_shares[3] # invasive only
  vary_lb_cf[i,7] = out_shares[4] # both tests
  vary_lb_cf[i,8] = out_shares[12] # terminated pregnancy without chrom abn
  vary_lb_cf[i,9] = out_shares[7] # live birth w CA & a >= c
  vary_lb_cf[i,10] = out_shares[17] - surplus_norm # consumer surplus
  
  rm(set, dec, out, out_shares)
}
print(vary_lb_cf)

write.csv(vary_lb_cf, file=paste0(RESULTS, "/vary_lb_cf.csv"), na=".", row.names=TRUE)

#### CA/preg share by bin ###
analysis_merge <- analysis %>%
  dplyr::select(pregnancy,bin_number)
all <- left_join(all, analysis_merge, by = "pregnancy")

all_ca <- all %>%
  filter(sim_positive == 1) %>% 
  dplyr::select(pregnancy, sim_wgt, p_i, fetus_risk, bin_number)

write.dta(all_ca, file=paste0(TEMP, "/all_ca.dta"))

print("End time")
end_time <- Sys.time()

print(end_time)
time_diff <- end_time - start_time
print(time_diff)

#### EXIT ####

sink()
#rm(list = ls())